//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Octant.hh
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Octant class definition.
 * \note   Copyright (C) 2011 Jeremy Roberts.  Derived from Denovo.
 */
//---------------------------------------------------------------------------//
// $Rev:: 133                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-08 19:31:48 +0000 (Thu, 08 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#ifndef OCTANT_HH
#define OCTANT_HH

#include <vector>

#include "angle/Quadrature.hh"
#include "angle/AngleTraits.hh"
#include "tran/Traits.hh"
#include "../utilities/Vector_Lite.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class Octant
 * \brief Store angles and weights by octant.
 *
 * The name "octant" is somewhat of a misnomer for application within
 * slabtran.  In 2-d and 3-d, it makes sense to speak of 4 and 8 octants as
 * comprising angular phase space.  However, in 1-d, we solve only in a the
 * polar angle \f$ \theta \f$, and there is really only a positive and
 * negative angular half-space; nonetheless, these half-spaces are represented
 * here as "octants".
 *
 * Construct an octant set of level-symmetric quadrature points on an octant.
 * The octant is defined by a triple \c (x-dir, \c y-dir, \c z-dir) coordinate
 * such that:
 * \f[
   \begin{array}{lll}
   \mathrm{octant}(+1, +1, +1) &\equiv &(+x, +y, +z)\:,\\
   \mathrm{octant}(-1, +1, +1) &\equiv &(-x, +y, +z)\:,\\
   \mathrm{octant}(+1, -1, +1) &\equiv &(+x, -y, +z)\:,\\
   \mathrm{octant}(-1, -1, +1) &\equiv &(-x, -y, +z)\:,\\
   \mathrm{octant}(+1, +1, -1) &\equiv &(+x, +y, -z)\:,\\
   \mathrm{octant}(-1, +1, -1) &\equiv &(-x, +y, -z)\:,\\
   \mathrm{octant}(+1, -1, -1) &\equiv &(+x, -y, -z)\:,\\
   \mathrm{octant}(-1, -1, -1) &\equiv &(-x, -y, -z)\:.
   \end{array}
 * \f]
 * These are defined in an Octant_Designator object.
 *
 * The number of angles in an octant is \f$ n = N(N+2) / 8\f$ for Level
 * Symmetric quadrature.  The quadrature set passed in must be a
 * Level-Symmetric set.  The quadrature set must be normalized to \f$4\pi\f$
 * over all weights.  Obviously, a single octant will not be normalized to
 * \f$ 4\pi\f$.
 *
 * The standard way to loop through angles and weights in an octant is:
 * \code
 *   Octant::angle_iterator angles = octant.angle_begin();
 *   Octant::weight_iterator wts   = octant.wt_begin();
 *   for (; angles != octant.angle_end(); angles++, wts++)
 *   {
 *       double wt           = *wts;
 *       Octant::Omega omega = *angles;
 *
 *       // do something ...
 *   }
 * \endcode
 *
 * The octants can also be defined to work in an adjoint mode.  In this mode,
 * the angles in each octant have opposite sign.  In forward mode the
 * following relationship holds:
 * \f[
   \begin{array}{lll}
   \mathrm{octant}(+1, +1, +1) &\rightarrow &(+\mu, +\eta, +\xi)\:,\\
   \mathrm{octant}(-1, +1, +1) &\rightarrow &(-\mu, +\eta, +\xi)\:,\\
   \mathrm{octant}(+1, -1, +1) &\rightarrow &(+\mu, -\eta, +\xi)\:,\\
   \mathrm{octant}(-1, -1, +1) &\rightarrow &(-\mu, -\eta, +\xi)\:,\\
   \mathrm{octant}(+1, +1, -1) &\rightarrow &(+\mu, +\eta, -\xi)\:,\\
   \mathrm{octant}(-1, +1, -1) &\rightarrow &(-\mu, +\eta, -\xi)\:,\\
   \mathrm{octant}(+1, -1, -1) &\rightarrow &(+\mu, -\eta, -\xi)\:,\\
   \mathrm{octant}(-1, -1, -1) &\rightarrow &(-\mu, -\eta, -\xi)\:,
   \end{array}
 * \f]
 * whereas in adjoint mode we have
 * \f[
   \begin{array}{lll}
   \mathrm{octant}(+1, +1, +1) &\rightarrow &(-\mu, -\eta, -\xi)\:,\\
   \mathrm{octant}(-1, +1, +1) &\rightarrow &(+\mu, -\eta, -\xi)\:,\\
   \mathrm{octant}(+1, -1, +1) &\rightarrow &(-\mu, +\eta, -\xi)\:,\\
   \mathrm{octant}(-1, -1, +1) &\rightarrow &(+\mu, +\eta, -\xi)\:,\\
   \mathrm{octant}(+1, +1, -1) &\rightarrow &(-\mu, -\eta, +\xi)\:,\\
   \mathrm{octant}(-1, +1, -1) &\rightarrow &(+\mu, -\eta, +\xi)\:,\\
   \mathrm{octant}(+1, -1, -1) &\rightarrow &(-\mu, +\eta, +\xi)\:,\\
   \mathrm{octant}(-1, -1, -1) &\rightarrow &(+\mu, +\eta, +\xi)\:.
   \end{array}
 * \f]
 * This allows sweeps in adjoint solves to go in the negative direction of
 * particles travel because adjoint solvers require calculating the incoming
 * angular flux given an outgoing flux.  This is opposite of forward solves in
 * which the outgoing angular fluxes are solved knowing the incoming values.
 *
 * The angles are ordered within each octant such that the reflected octant
 * has the same ordering as the incident octant at a face.  For example, the
 * reflected octant of \f$(+1, +1, +1)\f$ at the \f$i\f$-face, \f$(1, 0,
 * 0)\f$, is \f$(-1, 1, 1)\f$.  The angles have the same order within each
 * octant; thus, the \c m=0 angle in \f$(-1, 1, 1)\f$ is the reflected angle
 * of \c m=0 in \f$(1, 1, 1)\f$.  The angle of reflection is defined by
 * \f[
   \hat{\Omega}_{m'} = \hat{\Omega}_m - 2\hat{\mathbf{n}}_j
   (\hat{\Omega}_m\cdot\hat{\mathbf{n}}_j)\:.
 * \f]
 */
/*! 
 * \example angle/test/tstOctant.cc
 *
 * Test of Octant using an \f$S_8\f$ Level-Symmetric quadrature set.
 */
//===========================================================================//

template <class D>
class Octant
{

public:
  //@{
  //! Useful typedefs for Octant-dependent classes.
  typedef util::Vector_Lite<int, 3> Octant_Designator;
  typedef typename OmegaTraits<D>::Omega Omega;
  typedef std::vector<Omega> Angles;
  typedef std::vector<double> Weights;
  typedef std::pair<Omega, double> Angle_Weight;
  //@}

  //! Iterator for angles.
  typedef typename Angles::const_iterator angle_iterator;

  //! Iterator for weights.
  typedef Weights::const_iterator weight_iterator;

  // Constructor.
  Octant(const Octant_Designator &octant,
         const Quadrature<D> &quadrature,
         bool do_adjoint = false);

  // >>> INTERFACE

  //! Number of angles in the octant.
  int number_angles() const
  {
    return d_number_angles;
  }

  //! Get the octant designator.
  const Octant_Designator& get_octant() const
  {
    return d_octant;
  }

  //! Get the octant direction in a specified dimension.
  int direction(int ijk) const
  {
    return d_octant[ijk];
  }

  //! Get angles in this octant by subscript.
  const Omega& angles(int m) const
  {
    Require(v(m));
    return d_angles[m];
  }

  //! Get weights in this octant by subscript.
  double weight(int m) const
  {
    Require(v(m));
    return d_weights[m];
  }

  //! Return an iterator to the first angle in the octant.
  angle_iterator angle_begin() const
  {
    return d_angles.begin();
  }

  //! Return an iterator to the last angle in the octant.
  angle_iterator angle_end() const
  {
    return d_angles.end();
  }

  //! Return an iterator to the first weight.
  weight_iterator wt_begin() const
  {
    return d_weights.begin();
  }

  //! Return an iterator to the last weight.
  weight_iterator wt_end() const
  {
    return d_weights.end();
  }

  //! Query if this an adjoint octant.
  bool adjoint() const
  {
    return d_adjoint;
  }

  private:

  /// \name Data
  //\{

  //! Dimension
  int d_dimension;

  //! Number of angles per octant.
  int d_number_angles;

  //! Angles in this octant.
  Angles d_angles;

  //! Weights for each angle.
  Weights d_weights;

  //! Octant that these angles are defined in.
  Octant_Designator d_octant;

  //! Adjoint octants.
  bool d_adjoint;

  //\}

  /// \name Implementation
  //\{

  /// Validate an angle subscript.
  bool v(int m) const
  {
    return (m >= 0 && m < d_number_angles);
  }

  /// Load quadrature angles.
  void load_quadrature_angles(const Quadrature<D> &quadrature, double adj);

  //\}
};

} // end namespace slabtran

#endif // OCTANT_HH

//---------------------------------------------------------------------------//
//              end of Octant.hh
//---------------------------------------------------------------------------//
